1 Objetivo

Realizar a modelagem geoestatística da profundidade do solum área de produção silvicultural

2 Caracterização da área de estudo

Os dados que utilizei para este trabalho são provenientes de um povoamento florestal de 108 ha de Pinus taeda L. A área pertence ao município de Campo Belo do Sul, região serrana do Estado de Santa Catarina, Brasil. O clima é do tipo Cfb, mesotérmico, subtropical úmido e com precipitação média de 1.647 mm com chuvas bem distribuídas no ano. A geologia regional é constituída por uma sequência vulcânica de rochas ácidas da Formação Serra Geral, com predomínio de riodacito.

A área possui Neossolos Litólicos e Neossolos Regolíticos, Cambissolos Háplicos e Cambissolos Húmicos, Latossolos Vermelhos e Gleissolos Melânicos (Figura 1c). Conforme as tendências observadas no campo, em locais de relevo plano ou suavemente ondulado com boa drenagem estão solos profundos com sequência de horizontes A-Bw (Latossolos), em condições de má drenagem, ocorrem solos com a sequência de horizonte A-Cg (Gleissolos). Nos relevos ondulado ou fortemente ondulado, predominam solos rasos com a sequência de horizontes A ou A-Bi (Neossolos e Cambissolos).

Localização da área de estudo no município de Campo Belo do Sul, Estado de Santa Catarina (SC), Brasil (a) e área ampliada com modelo digital de elevação e distribuição dos pontos amostrais (b). Relação da classe de solo e altura das árvores de uma Topossequência típica da área de estudo (c).

Figure 2.1: Localização da área de estudo no município de Campo Belo do Sul, Estado de Santa Catarina (SC), Brasil (a) e área ampliada com modelo digital de elevação e distribuição dos pontos amostrais (b). Relação da classe de solo e altura das árvores de uma Topossequência típica da área de estudo (c).

3 Amostragem

A coleta de dados foi realizada em 94 pontos (Figura 1 b) alocados pelo algoritmo Hipercubo Latino Condicionado (Conditioned Latin Hipercube Sampling – cLHS) através da função clhs do pacote clhs para R. Foram consideradas como variáveis condicionantes da amostragem a elevação, profundidade do vale, índice de umidade topográfico, nivel base da rede de drenagem e declividade as quais, juntas, explicaram aproximadamente 86% da variância topográfica, identificada através da Análise de Componentes Principais (ACP).

Em cada ponto amostral foi medida a profundidade do solum (PF). Consideramos solum a espessura máxima do solo onde as raízes podem se desenvolver sem impedimentos físicos para penetração livre das raízes. Os fatores limitantes considerados foram o lençol freático elevado e o contato com rocha consolidada (contato lítico) com ou sem fissuras.

Cada ponto de amostragem representa uma área de aproximadamente 100 m2 no campo. Apesar do suporte amostral ser areal, considerei-os como sendo em ponto.

Os dados foram armazenados no objeto pontos definido como SpatialPolygonsDataFrame com as coordenadas projetadas em WGS84.

pontos <- read.csv('../data/GateadosDados.csv', dec = ".", sep= ";", stringsAsFactors = FALSE)
pontos$PFd <- pontos$PF/10
sp::coordinates(pontos) <- c('X' , 'Y')
wgs84utm22s <- sp::CRS('+proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
sp::proj4string(pontos) <- wgs84utm22s

A localização dos pontos no espaço pode ser observada a seguir:

mapview(pontos, zcol = "PFd")

4 Modelagem geoestatística da profundidade do Solum

4.1 Dados de profundidade do solum

Como uma análise preliminar avaliei a distribuição de frequências dos dados de profundidade do solum (Figura 2). Os dados apresentam uma distribuição bimodal, com valores entre 1 e 10 e média de 5,86 dm.
 Histograma de frequências da profundidade do solum

Figure 4.1: Histograma de frequências da profundidade do solum

summary(pontos$PFd)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   2.000   6.000   5.858  10.000  10.000

4.2 Modelo linear misto de variação espacial

Para modelagem geoestatística dos dados foi utilizado o modelo linear misto de variação espacial denotado por

\[Y(\boldsymbol{s}_i) = Z(\boldsymbol{s}_i) + \varepsilon(\boldsymbol{s}_i) = \boldsymbol{x}(\boldsymbol{s}_i)^\text{T}\boldsymbol{\beta} + B(\boldsymbol{s}_i) + \varepsilon(\boldsymbol{s}_i)\]

Para utilizar esse modelo foi necessário supor que os dados são uma realização de um campo aleatório \(Y(\boldsymbol{s}_i)\) com distribuição normal que podem ser descritos como a combinação aditiva de efeitos fixos, efeitos estocásticos e erro aleatório independente.

\(Z(\boldsymbol{s}_i)\) ou sinal possui dois componentes. O primeiro (efeito fixo) \(\boldsymbol{x}(\boldsymbol{s}_i)^\text{T}\boldsymbol{\beta}\) representa os efeitos de origem desterminística, que relaciona a dependência entre a variável e as covariáveis.

O segundo componente do sinal (efeito aleatório), \(B(\boldsymbol{s}_i)\), um campo aleatório Gaussiano estacionário não-observável, descrito por sua função de média e função de covariância.

\(\varepsilon(\boldsymbol{s}_i)\) é o erro (ou ruído), descrito por uma distribuição Gaussina de probabilidade, cujo parâmetro desconhecido de escala é \(\tau\).

4.2.1 Efeito fixo

Para selecionar os efeitos fixos do modelo, verifiquei a relação linear entre as variáveis topográficas e de produtividade e a profundidade do solum.

4.2.1.1 Covariáveis topográficas

As covariáveis topográfricas foram derivadas do modelo digital de elevação disponibilizado pelo Governo do Estado de SC - Secretaria de Estado do Desenvolvimento Econômico Sustentável, proveniente do Levantamento Aerofotogramétrico em 2010. Os dados, disponibilizados com resolução de 1 metro, foram reamostrados para resolução espacial de 10 metros utilizando a ferramenta reamostragem no software SAGA GIS.

Os planos de informação utilizados foram importados pela ferramenta raster::rastere armazenados em .tif nos objetos RasterLayer (Figura 3).

DECLI <- raster::raster("../data/Covars/DECLI.tif")
ELEV <- raster::raster("../data/Covars/ELEV.tif")
VD <- raster::raster("../data/Covars/VD.tif")
TWI <- raster::raster("../data/Covars/TWI.tif")
 Covariáveis topográficas utilizadas como preditoras no modelo floresta aleatória

Figure 4.2: Covariáveis topográficas utilizadas como preditoras no modelo floresta aleatória

A partir da função raster::extract extraí os valores de cada objeto na localização de cada observação contida no objeto espacial pontos.

pontos$DECLI <- raster::extract(DECLI, pontos)
pontos$ELEV <- raster::extract(ELEV, pontos)
pontos$VD <- raster::extract(VD, pontos)
pontos$TWI <- raster::extract(TWI, pontos)

4.2.1.2 Covariável de produção

Além dos atributos de terreno, foi utilizado um índice de produtividade disponibilizado pela empresa.

Esses índices são utilizados para ordenamento da produção e provém da amostragem de parcelas fixas. A área possui 12 parcelas fixas de inventário contínuo (PIC) de 500 metros quadrados cada que são utilizados na estimativa da produtividade local. Cada PIC é classificada em função da média de altura das 100 árvores de maior perímetro basal da parcela, expressa como altura dominante (Hdom). Em função da Hdom e seus incrementos anuais é estabelecido o índices de sítio (Sitio). Esse índice então é então atribuido à todo o talhão.

As parcelas da área são classificadas em 4 níveis, em que 1 corresponde a melhor e 4 a pior produtividade (Figura 4).

Utilizei a função raster::shapefile para carregar o polígono com informações das sítio - armazenado no objeto ProdutividadePistola. A função sp::spTransform foi usada para projetar as coordenadas original no plano cartesiano (UTM) e a função raster::extract para extraír os valores de cada objeto raster nas localizações de cada observação contidas no objeto espacial pontos.

Sitio <- raster::raster("../data/Covars/Sitio.tif")
Sitio <- as.factor(Sitio)
sp::proj4string(Sitio) <- wgs84utm22s
pontos$Sitio <- raster::extract(Sitio, pontos) %>%  as.factor()
pontos <- na.omit(pontos)
ProdutividadePistola <- 
 raster::shapefile('../data/Produtividade/ProdutividadePistola.shp') %>%
 sp::spTransform(wgs84utm22s)
ProdutividadePistola$Sitio <- as.factor(ProdutividadePistola$Sitio)

sp::spplot(
  ProdutividadePistola, scales = list(draw = T),
  main = 'Classificação dos sítios e  localização dos pontos') +
  lattice::xyplot(Y ~ X, data = as.data.frame(pontos@coords),
                  pch = 20, col = 'red', lwd = 2, cex = 1.5) %>%
latticeExtra::as.layer()
 Classificação dos sítios e  localização dos pontos na área de estudo

Figure 4.3: Classificação dos sítios e localização dos pontos na área de estudo

4.2.1.3 Avaliação da relação linear entre as variáveis

Considerando que apenas a covariável DECLI apresenta uma correlação significativa com a profundidade do solum, considerei somente a declividade como variável de efeito fixo, não-constantes, na área.

lm(PFd ~ DECLI + ELEV + VD + Sitio, data = pontos@data) %>% summary()
## 
## Call:
## lm(formula = PFd ~ DECLI + ELEV + VD + Sitio, data = pontos@data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.235 -2.569 -0.349  2.922  5.235 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  1.249e+01  4.116e+01   0.303  0.76223   
## DECLI       -1.517e+01  5.341e+00  -2.840  0.00557 **
## ELEV        -4.219e-03  4.458e-02  -0.095  0.92481   
## VD          -8.306e-04  5.096e-02  -0.016  0.98703   
## Sitio2       1.059e-01  1.089e+00   0.097  0.92270   
## Sitio3      -2.024e+00  1.105e+00  -1.832  0.07017 . 
## Sitio4      -1.257e+00  9.954e-01  -1.263  0.20998   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.381 on 91 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.1666, Adjusted R-squared:  0.1116 
## F-statistic: 3.031 on 6 and 91 DF,  p-value: 0.009525

4.3 Variograma amostral

O variograma amostral (Figura 7) foi computado através da função georob::sample.variogram. O estimador para semivariãncia foi Metheron (método dos momentos). Para a obtenção dos parâmetros utilizei um corte de 66% da distância máxima entre os pontos, armazenada no objeto distmax, excluindo os pares de longo alcance.

distmax <-dist(pontos@coords) %>% max() /3

limites <- seq(0, distmax, length.out = 20)
vario<- georob::sample.variogram(PFd ~ DECLI,
    data= pontos, locations = ~ X + Y, lag.dist.def = limites, estimator = "matheron") %>%
plot(ylab = 'Seminvariância', xlab = 'Distância de separação (m)', annotate.npairs = TRUE, main = "Semivariograma")
 Variograma amostral

Figure 4.4: Variograma amostral

4.4 Ajuste da função que descreve o varigrama amostral

O método de ajuste empregado no variograma amostral foi o de quadrados mínimos não-lineares ponderados, com ajuste de um modelo exponencial, com ponderação definida conforme o método de “Cressie”. O processo de estimativa dos parâmetros do modelo exponencial do variograma foi conduzido via otimização usando a função stats::optim(method = "BFGS"). A função resultante desse processo foi armazenada no objeto vario_fit.

vario_fit <- 
  georob::fit.variogram.model(
  vario, variogram.model = 'RMexp', param = c(variance = 10, nugget = 5, scale = 70), weighting.method = "cressie", method = "BFGS")
## 
## 
##                             variance       snugget        nugget         scale
##   Variogram parameters 7.666473e-122  0.000000e+00  1.458980e-65  9.136816e+14
summary(vario_fit)
## 
## Call:georob::fit.variogram.model(sv = vario, variogram.model = "RMexp", 
##     param = c(variance = 10, nugget = 5, scale = 70), weighting.method = "cressie", 
##     method = "BFGS")
## 
## Convergence in 97 function and 34 Jacobian/gradient evaluations
## 
## Residual Sum of Squares: 58.86111 
## 
## Residuals (epsilon):
##     Min      1Q  Median      3Q     Max 
## -4.0586 -1.8355 -0.6017  1.5709  2.8500 
## 
## Variogram:  RMexp 
##                Estimate    Lower Upper
## variance       12.43754  9.10081 17.00
## snugget(fixed)  0.00000       NA    NA
## nugget          1.28807  0.06083 27.28
## scale          72.24013 52.67343 99.08

O ajuste do modelo exponencial ao variograma amostral é mostrado na Figura 8. A curva ajustada passa próximo ao centro de massa dos vinte pontos do variograma amostral.

plot(vario, xlab = 'Distância de separação (m)', ylab = 'Semivariância', annotate.npairs = TRUE)
lines(vario_fit, col = "red", lty = 'dashed')
 Variograma amostral (em preto) e função exponencial a ele ajustado (vermelho)

Figure 4.5: Variograma amostral (em preto) e função exponencial a ele ajustado (vermelho)

Considerando que durante a obteção dos dados de campo as informações de profundidade do solum os valores foram obtidos em decímetros, considerei que a variância do erro de medida é igual a 0,5 dm devido ao arredondamento dos valores.

Assim, foi possível discretizar a variância não explicada em erro de medida corresponde 0,25 do parâmetro nugget. A variância restante foi atribuída a variação espacial não auto-correlacionada espacialmente - não capturada pelo plano amostral snugget. A função resultante desse processo foi armazenado no objeto vario_fit_error.

nugget <- 0.25

vario_fit_error <- georob::georob(
   PFd ~ DECLI, pontos, locations = ~ X + Y, variogram.model = 'RMexp', 
 param = c(variance = vario_fit$variogram.object[[1]]$param[['variance']], 
           nugget = nugget,
           snugget = vario_fit$variogram.object[[1]]$param[['nugget']] - nugget,
           scale = vario_fit$variogram.object[[1]]$param[['scale']]),
 fit.param = georob::default.fit.param(nugget = FALSE, snugget = TRUE),
 tuning.psi = 1000, control = georob::control.georob(initial.fixef = 'lm'))
summary(vario_fit_error)
## 
## Call:georob::georob(formula = PFd ~ DECLI, data = pontos, locations = ~X + 
##     Y, variogram.model = "RMexp", param = c(variance = vario_fit$variogram.object[[1]]$param[["variance"]], 
##     nugget = nugget, snugget = vario_fit$variogram.object[[1]]$param[["nugget"]] - 
##         nugget, scale = vario_fit$variogram.object[[1]]$param[["scale"]]), 
##     fit.param = georob::default.fit.param(nugget = FALSE, snugget = TRUE), 
##     tuning.psi = 1000, control = georob::control.georob(initial.fixef = "lm"))
## 
## Tuning constant:  1000 
## 
## Convergence in 8 function and 6 Jacobian/gradient evaluations
## 
## Estimating equations (gradient)
## 
##                                   xi         scale
##   Gradient           :  1.069608e-03 -8.586508e-04
## 
## Maximized restricted log-likelihood: -263.634 
## 
## Predicted latent variable (B):
##     Min      1Q  Median      3Q     Max 
## -5.8300 -2.9564 -0.1296  3.6171  5.4093 
## 
## Residuals (epsilon):
##        Min         1Q     Median         3Q        Max 
## -0.1847714 -0.0623897  0.0009239  0.0657684  0.1921513 
## 
## Standardized residuals:
##      Min       1Q   Median       3Q      Max 
## -2.29890 -0.79556  0.01281  0.88223  1.98182 
## 
## 
## Gaussian REML estimates
## 
## Variogram:  RMexp 
##                Estimate  Lower Upper
## variance         10.696  7.607 15.04
## snugget(fixed)    1.038     NA    NA
## nugget(fixed)     0.250     NA    NA
## scale            48.614 24.323 97.16
## 
## 
## Fixed effects coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.2931     0.7609   9.585 8.05e-16 ***
## DECLI       -12.1389     5.0173  -2.419   0.0174 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error (sqrt(nugget)): 0.5 
## 
## Robustness weights: 
##  All 102 weights are ~= 1.

A comparação entre as funções vario_fit_error e vario_fit está na figura 9.

plot(vario)
lines(vario_fit, col = "blue")
lines(vario_fit_error, col = "red")
 Variograma amostral (em preto) do modelo linear da profundidade do solum e a função exponencial a ele ajustada (azul) e a função ajustada com erro de medida fixo

Figure 4.6: Variograma amostral (em preto) do modelo linear da profundidade do solum e a função exponencial a ele ajustada (azul) e a função ajustada com erro de medida fixo

A função vario_fit_error foi utilizada para a predição espacial. Para isso foi criado um grid de predição abrangendo a área de estudo. A partir da função raster::extract foram extraídos os valores do objeto DECLI e armazenados no novo conjunto de dados grid.

grid <- sp::spsample(ProdutividadePistola, 10000, type = 'regular')
grid$DECLI <- raster::extract(DECLI, grid)

grid<- 
  sp::SpatialPointsDataFrame(
    coords = grid@coords, 
    data = data.frame(grid),
    proj4string = grid@proj4string)
colnames(grid@coords) <- colnames(pontos@coords)

4.5 Predição espacial da resposta do modelo linear misto

Utilizei a função predict para aplicar a função vario_fit_erroraos pontos de toda área armazenados no objeto grid

pred_ponto <- predict(
  vario_fit_error, newdata = grid, type= "response", signif = 0.95,
  control = georob::control.predict.georob(extended.output = TRUE))
sp::gridded(pred_ponto) <- TRUE
spplot(pred_ponto)
Mapas de predição - saída extendida do georob

Figure 4.7: Mapas de predição - saída extendida do georob

trend é a média espacial dos efeitos deterministicos (declividade), suavisados.

A krigagem tenta minimizar o valor esperado do erro entre o valor predito e o medido. Ainda assim o intervalo de predição ainda é grande

A raiz do SE é erro padrão da krigagem. O quadrado do SE é a variância dos erros de predição.

at <- pred_ponto@data[, c("pred", "lower", "upper")] %>% range(na.rm = TRUE)
at <- seq(at[1], at[2], length.out = 20)

spplot(pred_ponto, zcol="lower", at=at, main="Limite inferior")

spplot(pred_ponto, zcol="pred", at=at, main="Realização mais provável")

spplot(pred_ponto, zcol="upper", at=at, main="Limite superior")

4.5.1 Avaliação da qualidade da predição

Para a validação dos resultados utilizei o método de validação cruzada. Este método é realizado em etapas de validação com a partição do conjunto de dados aleatóriamente em subconjuntos (\(k\)). Em cada etapa de validação um dos conjuntos é utilizado para validar o modelo calibrado com \(k-1\) subconjuntos. Esse procedimento é repetido até que todos subconjuntos \(k\) tenham sido utilizados como conjunto de validação do modelo. A partir das predições realizadas com cada subconjunto \(k\) são calculados os erros para avaliar a qualidade das predições.

Para isso utilizei a função cv::georob. Defini o número o número \(k\) de subconjuntos em que o conjunto de dados pontos foi particionado definindo o parâmetro nset = 101. Assim, em cada etapa da validação 101 pontos foram utilizados para a validação e 1 para a calibração do modelo linear misto.

validacao <- georob::cv(vario_fit_error, nset = 101)
summary(validacao)
## 
## Statistics of cross-validation prediction errors
##       me      mede      rmse      made       qne      msse    medsse  
## -0.04291   0.03903   3.29989   4.27377   3.35070   1.08072   0.76194  
##     crps  
##  1.92458
1 - sum((validacao$pred$data - validacao$pred$pred)^2) / sum((validacao$pred$data - mean(validacao$pred$data))^2)
## [1] 0.1323907
plot(validacao)

Como os dados utilizados para a validação foram os mesmo utilizados para validar o modelo, e, que as amostras foram obtidas intencionalmente, é provavel que esta avaliação seja otimista.